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1  Introduction 

In  this  project  we  have  developed  some  of  the  necessary  tools  and  soft¬ 
ware  infrastructure  to  enable  multiscale  approaches  for  electronic  transport 
and  heat  dissipation  in  materials  and  devices.  The  multiscale-multiphysics 
concept  requires  exchange  of  information  between  the  microscopic  and  the 
macroscopic  levels  of  descriptions,  allowing  simulations  of  the  physical  be¬ 
havior  of  a  material  including  nanoscale  details.  This  computational  paradigm 
allows  to  bridge  over  several  orders  of  magnitudes  of  scale-lengths  and  time- 
scales,  transferring  information  between  the  micro  and  the  macro  world  or 
vice-versa.  The  multiscale  and  multi  physics  approach  can  be  applied  to  the 
most  diverse  physical  problems  and  disciplines  from  biochemistry  to  ma¬ 
terials  science.  Examples  of  current  working  applications  range  from  the 
chemical  behavior  of  molecular  reactions  to  protein  and  enzyme  functions 
or  ion  pumps,  the  analysis  of  structural  defects  and  crack  formations  to 
problems  of  surface  adhesion  and  cathalisis.  Successful  applications  have 
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also  been  reported  in  the  study  of  photoexcitations,  exciton  dissociations 
and  charge  transfer.  In  most  cases  the  multiscale  approach  consists  in  com¬ 
bining  quantum  mechanical  calculations  (QM)  with  faster  semi-empirical  or 
empirical  interatomic  forces  (molecular  mechanics).  In  the  last  decade  sev¬ 
eral  attempts  have  been  made  to  couple  interatomic  forces,  obtained  with 
semiempirical  or  empirical  potentials,  with  macroscale  simulations,  usually 
performed  using  finite-elements  schemes,  describing  the  materials  with  av¬ 
erage  local  parameters  such  as  elastic  properties,  electron/ion  mobilities, 
energy  bands,  etc.  This  is  the  last  chain  of  a  whole  hierarchy  of  meth¬ 
ods,  but  hardly  easy  to  accomplish.  Coupling  macroscopic  with  atomistic 
models  poses  several  difficulties  due  to  quite  heterogeneous  frameworks  and 
formalisms  involved,  since  the  two  approaches  have  been  historically  devel¬ 
oped  by  different  communities  to  target  substantially  different  problems. 

Multiscale  approaches  can  be  applied  both  to  obtain  macroscopic  pa¬ 
rameters  from  microscopic  models  ( overlap  method)  or  different  scales  can 
be  sewed  together  at  some  boundary  ( bridge  method)  in  order  to  treat  criti¬ 
cal  subregions  with  appropriate  models.  Both  shemes  are  valid  alternatives, 
depending  on  the  problem  to  be  solved,  but  bridge  methods  require  mi¬ 
cro  and  macro  scales  to  tightly  couple  and  exchange  information  in  both 
directions  posing  additional  interfacing  challenges.  Practical  examples  of 
this  methodology  are  the  multiscale  models  developed  to  compute  strain 
relaxation  around  a  crystal  defect  or  to  nanoelectromechanical  cantilevers. 

The  processes  of  interest  in  this  project  are  electron  and  heat  transport 
in  materials  and  devices,  where  the  multiscale  concept  has  not  been  fully 
developed  yet.  Recent  technological  advances  have  progressively  scaled  the 
active  region  of  electronic  and  optical  devices  down  to  scales  of  few  tents  of 
nanometer  or  lower.  However  the  nano- world  is  always  addressed  via  macro¬ 
scopic  contacts  that  carry  the  current  information.  Heat  sinks  are  also  of 
macroscopic  size  and  heat  managment  has  becoming  a  pressing  problem  in 
modern  nanoelectronics.  It  is  therefore  of  crucial  importance  for  the  devel¬ 
opment  of  the  field  to  couple  micro  and  macro  scales  in  order  to  address 
electro-thermal  properties.  Possible  applications  of  this  multiscale  coupling 
range  from  electrical  simulations  of  nanodevices  and  their  thermal  manag¬ 
ment.  The  model  will  also  be  applicable,  with  suitable  adaptations,  to  study 
heat  dissipation  in  materials  where  heat  is  locally  produced  with  mechanisms 
other  than  joule  heating,  for  instance  as  a  consequence  of  light  absorption. 

The  whole  point  of  this  research  is  to  build  a  suitable  environement  in 
which  electronic  and  heat  transport  simulations  (performed  at  FEM  level) 
can  be  coupled  with  microscopic  models  for  heat  generation  and  transport. 
The  heat  generation  rate  at  the  micro  scale  can  be  accounted  by  imposing 
corresponding  boundary  conditions  at  the  micro-macro  boundaries. 

Among  the  several  available  software  implementing  FE  calculations,  we 
have  chosen  the  LibMesh  library,  which  is  a  well  established  and  lively  open- 
source  (LGPL)  project.  Both  aspects  are  quite  important  when  developing 
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new  derivative  works,  since  the  feedback  from  library  authors  for  clarifi¬ 
cations  and  bug  fixes  may  be  necessary  in  many  cases.  LibMesh  contains 
interfaces  to  the  PETSc  and  SLEPc  numerical  solvers  and  the  METIS  par¬ 
titioning  library,  giving  a  complete  environment  for  development  of  general 
solutions  of  partial-differential  equations  based  on  the  FEM  approach. 


2  Heat  dissipation  models 

2.1  The  Fourier  transport  equations 

Heat  balance  problems  at  the  macroscopic  scale  are  governed  by  the  Fourier’s 
law,  which  together  with  the  continuity  equation  reads  as: 


V  •  (— kVT)  =  Hs, 


(1) 


where  ft  is  the  lattice  thermal  conductivity  tensor,  T  the  temperature  and 
Hs  is  the  total  heat  source.  The  heat  source  may  be  given  by  any  dissipa¬ 
tive  processes  such  as  light  absorption  or  Joule’s  effect.  The  latter  can  be 
associated  to  dissipation  induced  by  carrier  transport  and  can  be  described 
at  macroscopic  scales  in  terms  of  a  drift-diffusion  model.  The  heat  source 
can  be  written  explicitly  in  terms  of  the  electron  and  hole  currents  and 
thermoelectric  effect, 


Hs  —  —V  •  [( PnT  +  </>n)Jn  +  (PpT  +  <t>p) Jp]  5 


(2) 


where  <fn  and  <fp  are  the  electro-chemical  potentials  and  Pn  and  Pp  are  the 
thermoelectric  power  for  electrons  and  holes  given  by 


and 


_  kb  /  5  e(j)n  +  Ec  -  e(p 
n~  q\  2  kbT 


p  _  h  (5  _e<Pp  + Ev -e<p 


q  \2  kbT 

The  expression  2  can  be  split,  as  shown  in  Table  1 


(3) 

(4) 


Expression 

Heat  source 

|Jn|'2 

&  n 

|JPI2 

qRsRH  (4>p  —  4>n  +  T  (Pp  —  Pn )) 

TJn  •  V Pn 
-TJp  •  VPp 

Electrons  Joule’s  effect 

Holes  Joule’s  effect 

Recombination  effect 

Electrons  Peltier- Thoms  on  effect 
Holes  Peltier- Thomson  effect 

Table  1:  Drift  diffusion  heat  sources 

The  equations  above  have  been  implemented  on  a  finite-element  (FEM) 
solver  based  on  the  LIBMESH  library.  This  solver  is  part  of  the  Tiber  CAD 


4 


simulator  under  development  in  the  ’Tor  Vergata’  group.  The  project  can 
be  visioned  at  the  mi  http://www.tibercad.org. 

As  a  working  example  we  report  a  2-dimensional  simulation  of  a  PN 
diode  surrounded  by  air.  The  finite  element  meshing  of  the  device  is  shown 
in  Fig.  1.  The  temperature  at  the  air  boundary  was  fixed  to  300  K  and 


Figure  1:  (Left)  Mesh  structure  and  (Right)  temperature  map  of  a  pn  diode 


we  have  included  a  thermal  surface  resistance  representing  heat  dissipation 
into  the  substrate.  The  temperature  map  reported  in  Figure  1  refers  to  a 
direct  bias  of  the  diode  and  shows  a  rise  of  temperature  near  the  pn  junction 
where  e/h  recombination  takes  place. 

Different  types  of  boundary  conditions  can  be  handled  already  within 
our  tool,  namely 

•  1.  Insulating  surface:  Jq  •  N  =  0 

•  2.  Ideal  conducting  surface:  T  =  Text 

•  3.  Thermal  resistive  surface:  Jq  •  N  =  GS(T  —  Text) 

In  the  above,  Jq  is  the  thermal  flux  and  Gs  is  the  surface  thermal  conduc¬ 
tivity. 

As  an  example  of  heat  generators  we  have  modeled  a  simple  array  of  dots 
acting  as  heat  sources.  The  Fourier  transport  model  is  solved  by  imposing  a 
fixed  temperature  at  the  boundary  and  a  realistic  source  of  heat  emanating 
from  each  dot.  The  results  are  shown  in  figure  2. 

2.2  Microscopic  theory  of  heat  transport:  grey  model 

As  the  characteristics  length  becomes  comparable  with  the  phonon  mean 
free  path  the  approximation  of  local  thermal  equilibrium  introduces  a  not 
negligible  error.  Several  models  based  on  the  Boltzmann  Transport  Equa¬ 
tion  (BTE)  have  already  been  proposed  in  the  last  decade  [1].  We  have 
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Figure  2:  (Left)  Array  of  heat  sources  (Right)  Temperature  map 


explored  applications  of  the  so-called  grey  model  [2],  which  relies  on  the 
following  assumptions:  (i)  collision  experienced  by  phonons  are  described 
via  a  relaxation  time,  which  does  not  depend  on  the  phonon  energy;  (ii)  All 
phonons  are  considered  to  have  the  same  group  velocity,  equal  to  the  first 
sound  velocity.  Neglecting  the  full  phonon  dispersion  is  a  severe  assumption 
that  will  be  dropped  in  future  developments.  With  these  assumptions  the 
BTE  for  phonons  is 

vgV  ■  ss  =  +  H  (5) 

T 

where  e  is  the  density  of  phonon  energy  propagating  along  the  s  direction, 
vg  is  the  sound  velocity,  eo  the  equilibrium  energy  density,  r  the  scattering 
relaxation  time  and  H  the  total  heat  source  per  solid  angle.  Energy  conser¬ 
vation  is  achieved  by  setting  a  global  condition  on  the  energy  flux  over  the 
whole  solid  angle.  As  equation  (5)  needs  to  be  solved  for  each  direction,  s, 
we  need  to  discretize  the  solid  angle.  The  thermal  flux  in  direction  i  can  be 
computed  using 

J\  =  i^Jn  v9£^dn-  (6) 

Finally  we  relate  the  local  lattice  temperature  with  the  density  of  heat  energy 
using  [3] 

T  =  T0  +  (e-eo)^ry  (7) 

where  C(Tq )  is  the  heat  capacity  at  the  temperature  Tq.  The  all  algorithm 
can  be  summarized  as  follows 


•  1.  Set  the  first  guess  as  the  equilibrium  energy  density,  eo  and  isotrop¬ 
ically  spread  it  in  all  directions. 
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•  2.  Solve  Eq.  (5)  for  each  direction. 

•  3.  Compute  the  new  equilibrium  energy. 


•  4.  Go  to  step  2  until  convergence  is  reached. 


Dirichlet  boundary  conditions  can  be  employed  by  fixing  the  generated 
phonon  energy  to  the  desired  value.  Thermal  insulating  boundary  condition 
can  be  applied  by  considering  both  specular  and  diffusive  interfaces.  The 
error  is  calculated  for  the  equilibrium  density.  It  is  possible  to  show  that 
whenever  the  mean  free  path  becomes  much  smaller  than  the  device  size, 
the  gray  model  approaches  the  diffusive  limit  and  Fourier  thermal  transport 
is  recovered. 

2.3  Towards  the  coupling  of  micro- macro  scales 

In  this  section  we  show  a  practical  coupling  of  heat  transport  models  at 
the  micro  (grey  model)  and  macro  (Fourier)  scales.  To  this  end  we  need  to 
define  a  general  embracing  strategy  between  microscopic  and  macroscopic 
domains  and  define  their  exchange  of  information.  We  need  to  distinguish 
between  two  cases  of  multiscale  couplings.  The  first  situation  is  when  the 
micro  and  macro  regions  spatially  coincide.  We  may  call  this  a  overlapping 
condition.  This  situation  could  be  useful  to  investigate  the  properties  of 
molecules  embedded  in  hosting  materials,  such  as  polymers  or  nanostruc- 
tured  media,  and  may  be  applied  especially  to  study  materials  coupled  to 
light-sources,  like  OLEDs,  lasers  or  energy  conversion  devices  such  as  or¬ 
ganic  solar  cells.  In  this  case  the  microscopic  calculations  will  couple  to 
the  classical  electro-thermal  equations  in  the  form  of  source  terms  that  will 
reflect  the  exciton  generation  recombination  and  dissociation  rates  and  will 
also  provide  power  released  per  unit  volume  under  the  form  of  heat.  In  the 
second  case,  that  we  may  call  bridging ,  the  microscopic  region  is  embedded 
within  the  classical  (macroscopic)  regions  and  the  two  are  spatially  sepa¬ 
rated.  In  this  case  the  two  regions  communicate  via  appropriate  continuity 
conditions  of  the  heat  flux  across  the  boundary  between  the  two  regions. 
The  latter  will  come  from  the  calculation  of  the  rate  of  energy  flowing  at  the 
interface.  The  models  in  the  two  regions  will  then  be  solved  self-consistently 
by  implementing  an  appropriate  feedback  loop.  To  many  respects  this  sit¬ 
uation  is  more  difficult  to  handle  than  the  first,  since  the  convergence  of 
the  macroscopic  FEM  equations  and  application  of  flux  boundaries  is  more 
delicate.  For  this  reason  the  bridge  case  has  been  here  analyzed  in  detail  in 
order  to  asses  the  feasibility  of  the  approach. 


7 


In  principle,  we  might  perform  the  BTE  over  the  whole  domain,  but 
the  discretization  over  the  solid  angle  would  strongly  increase  the  compu¬ 
tational  effort.  For  this  reason  we  compute  the  BTE  only  in  the  region 
where  it  is  needed  (micro  domain)  and  the  Fourier  model  is  applied  to  the 
rest  of  the  device  (macro  domain).  Both  models  can  communicate  through 
their  boundary  conditions  in  a  bridge  scheme.  The  Fourier  domain  fixes  the 
boundary  temperature  to  the  grey  model.  Conversely,  the  latter  fixes  the 
heat  flux  to  the  Fourier  model.  The  self-consistent  scheme  that  we  have 
developed  can  be  summarized  as  follows 


•  1.  Solve  the  Fourier  model  on  the  whole  domain  (micro  and  macro) 
which  fixes  the  boundary  temperature  on  the  micro  domain. 

•  2.  Solve  the  gray  model  on  the  micro  domain  which  imposes  the 
boundary  heat  flux  on  the  macro  domain. 

•  3.  Solve  the  Fourier  model  on  the  macro  domain. 

•  4.  Go  to  step  2  if  the  convergence  is  not  reached. 


In  order  to  speed-up  convergence  we  find  that  a  crucial  ingredient  is  the 
first  guess  of  the  equilibrium  energy  density.  A  good  guess  was  obtained 
by  computing  the  Fourier  heat  transport  over  the  whole  device,  including 
the  microscopic  region.  On  the  other  hand  we  find  that  whenever  the  grey 
domain  becomes  larger  than  the  phonon  mean  free  path,  the  number  of  steps 
required  for  convergence  significantly  increases. 

2.4  Application  to  a  GaN-QD  embedded  in  AlGaN 

The  method  is  applied  on  a  system  comprising  a  5  nm  GaN  cubic  quantum 
dot  embedded  in  a  nanocolumn  of  Alo.2GaN  of  100  nm  in  height.  The  two 
contacts  are  p  and  n  doped  l19  cm-3,  respectively.  The  QD  is  surrounded 
by  a  buffer  layer  of  intrinsic  Alo.2GaN  (see  Fig.  1).  We  first  compute  the 
I-V  characteristics  of  the  device  within  the  drift  diffusion  model,  showing  a 
diode  behaviour,  as  expected.  Joule  heat  is  generated  by  holes  and  electron 
transport,  as  shown  in  Figure  3  The  dominant  contribution  of  heat  due  to 
recombination  within  the  QD.  The  solid  angle  for  the  integration  of  the  grey 
model  equation  is  discretized  in  18  parts  (6  the  azimut  angle  and  3  the  ele¬ 
vation  angle).  Sounds  velocity  and  thermal  conductivities  of  GaN  and  AIN 
are  taken  from  the  literature  [4].  Thermal  relaxation  times  are  computed 
from  the  relation  3 k/  ( Cv2g )  in  order  to  be  consistent  with  the  diffusive  limit 


Figure  3:  (Left)  Structure  of  a  Nanocolum  with  a  QD.  (Right)  heat  sources  due  to 
n  and  p  joule  effects  and  recombination  heat 


[3].  The  sound  velocity  of  the  AlGaN  is  obtained  by  relying  on  the  virtual 
crystal  approximation  as  in  [4].  As  boundary  conditions  at  end  of  the  two 
contacts  we  set  the  temperature  of  300  K.  As  shown  in  Fig.  4  the  BTE  and 
the  Fourier  model  give  a  maximum  temperature  of  about  334  K  and  319 
K,  respectively,  revealing  that  that  phonons  distribution  is  far  away  from 
the  local  equilibrium.  Furthermore,  the  BTE  model  gives  a  strongly  peaked 
temperature  across  the  quantum  dot  region,  resulting  in  a  more  realistic 
temperature  profile.  Lets  now  consider  the  multiscale  approach.  The  grey 
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Figure  4:  (Left  panel)  Comparison  of  temperature  profiles  from  Fourier  and  gray 
model.  (Right  panel)  Multiscale  coupling  of  Fourier  and  grey  models. 

model  is  applied  to  the  micro  domain  which  includes  the  quantum  dot  and 
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the  buffer  regions,  whereas  the  contacts  belong  to  the  macro  domain,  treated 
with  the  Fourier  model.  The  two  models  are  solved  self-consistently  as  out¬ 
lined  above.  The  temperature  profile,  as  shown  in  the  right  panel  of  Fig.  4, 
matches  the  BTE  model  performed  over  the  whole  domain.  In  fact  one  may 
expect  at  first  that  the  multiscale  temperature  profile  should  matches  the 
Fourier  solution  in  the  regions  where  Fourier  is  solved.  However,  we  should 
bare  in  mind  that  the  multiscale  solution  knows  about  the  BTE  solved  in 
the  QD  region,  which  leads  to  a  different  effective  boundary  condition  for 
the  temperature  at  the  micro/macro  interface,  reflected  in  a  different  tem¬ 
perature  profile.  On  the  other  hand,  the  Fourier  model  and  the  gray  model 
correctly  approach  each  other  at  length-scales  larger  than  the  phonon  mean 
free  path.  We  finally  point  out  that  the  size  of  the  micro  region  is  a  crucial 
choice  that  depends  on  the  phonon  mean  free  path  and  the  lattice  thermal 
conductivity,  as  well  as  the  magnitude  of  the  heat  source.  Since  the  tem¬ 
perature  affects  the  energy  gap,  the  optical  recombination  spectrum  can  be 
influenced  by  the  heat  transport  model.  As  shown  in  Fig.  5  the  energy  peak 
undergoes  a  red  shift  of  about  10  meV  when  the  BTE  is  employed  in  place 
of  the  Fourier  model. 


Figure  5:  Optical  recombination  spectrum  in  the  QD  obtained  for  different  thermal 
models. 


2.5  Conclusions 

A  multiscale  scheme  for  thermal  modeling  of  nanostructures  has  been  devel¬ 
oped  and  applied  to  a  GaN  quantum  dot  embedded  in  a  AlGaN  nanocolumn. 
Temperature  results  of  Fourier,  BTE  and  multiscale  models  are  compared. 
We  found  out  that  the  Fourier  model  underestimates  the  temperature  across 
the  quantum  dot  and  justifies  the  need  of  a  BTE  based  model.  The  max¬ 
imum  temperature,  across  the  dot,  is  about  319  K  for  the  Fourier  model, 
against  334  K  of  the  BTE  model.  We  further  identified  a  convenient  choice  of 
the  initial  guess  of  equilibrium  energy  density,  which  strongly  speeds  up  the 
loop  convergence  of  the  gray  model.  Finally,  we  investigated  the  influence 
of  the  used  model  on  the  optical  spectrum. 
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3  Heat  dissipation  in  molecular  junctions 


The  first  step  of  development  consists  in  formulating  a  microscopic  theory 
of  heat  dissipation  due  to  joule  heating  in  nanoj unctions.  We  focus  our 
attention  to  molecular  bridges  as  for  instance  those  obtained  in  nanopores, 
break-j unctions  or  STM  configurations.  This  description  can  be  later  ex¬ 
tended  to  other  types  of  electronic  nanoj  unctions,  such  as  semiconducting 
quantum  dots  or  quantum  wells. 

At  the  microscopic  level,  the  release  of  energy  in  a  molecular  junction  is 
described  by  a  quantum  mechanical  approach,  fully  including  the  electron- 
phonon  coupling.  Molecular  vibrations  and  electron-vibron  couplings  are 
treated  in  the  usual  way  [5,  6,  7,  8],  after  decoupling  the  Hamiltonian  as 
a  superposition  of  independent  one-dimensional  oscillators  corresponding  to 
the  normal  modes  of  vibrations,  labeled  by  q.  The  harmonic  oscillators  can 
be  quantized  following  the  usual  prescriptions,  by  making  use  of  the  standard 
relationships  between  the  position  operator  and  the  creation/annihilation 
field  operators,  a+  and  aq.  Therefore,  the  ionic  vibrations  of  a  molecular  wire 
can  be  described  as  a  class  of  bosonic  particles,  here  referred  as  phonons .  The 
result  is  a  term  in  the  Hamiltonian  describing  the  electron-phonon  coupling, 


H-el—ph  ^  ^ 


7 11VC11  cv 


an  +  aQ 


(8) 


where  c+  and  cv  are,  respectively,  the  creation  and  annihilation  operators  of 
one  electron  in  the  atomic  basis-set  and  7^  are  the  electron-phonon  coupling 
matrices.  Mg  are  the  atomic  masses,  uq  the  mode  frequencies  and  e are  the 
normalized  atomic  displacements  for  each  mode.  The  non-orthogonality  of 
the  atomic  basis-set  is  reflected  by  the  presence  of  the  overlap  matrix, 
and  its  derivative  with  respect  to  the  ionic  positions,  Ka. 

While  crossing  the  system,  electrons  interact  with  the  molecular  ionic 
vibrations  from  which  they  can  be  inelastically  scattered.  The  formalism 
of  our  choice  is  the  non-equilibrium  Green’s  function  approach,  which  is 
very  well  suited  to  study  such  problems.  The  electron-phonon  interaction  is 
treated  within  perturbation  theory  of  the  non-equilibrium  Green’s  function 
formalism  and  the  current  through  the  junction  is  computed  using  the  Meir- 
Wingreen  formula  [9]  which  contains  an  elastic  and  an  inelastic  term.  The 
electron-phonon  self-energy  can  be  evaluated  with  diagrammatic  techniques 
such  as  the  self-consistent  Born  approximation  (SCBA). 


Spfc(>)  M  =  *  T  fAgG<(>)  (u  -  u/)j qD<>^  (u/) ,  (9) 

where  Dq  ^  are  the  correlation  functions  related  to  the  vibrational  modes. 
Explicit  forms  for  the  self-energy  can  be  found  when  the  vibrations  are 
assumed  Einstein  oscillators  with  delta-like  density  of  states.  Such  forms 
are  presented  in  [10,  11]. 
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3.1  Phonon  rate  equation 


When  a  bias  is  applied  and  a  current  starts  to  flow,  the  interaction  between 
electrons  and  vibrations  drives  the  population  of  phonons  out  of  equilibrium. 
That  is,  the  distribution  of  phonons  deviates  from  a  simple  Bose-Einstein 
distribution  due  to  phonon  emission  and  absorption  processes.  Moreover, 
the  presence  of  the  leads  allows  the  relaxation  of  phonons  from  the  device 
into  the  contacts.  In  order  to  describe  these  processes,  a  semi-classical  rate 
equation  has  been  implemented  [10], 

dE  =  Rq-Jq[Nq-nq{T0)],  (10) 

where  Nq  is  the  non-equilibrium  phonon  population,  nq(To )  is  the  equi¬ 
librium  phonon  distribution  at  temperature  To,  and  Rq  is  the  net  rate  of 
phonons  emitted.  This  rate  is  computed  from  the  emitted  power  in  each 
mode,  Pg,  which  can  be  computed  according  to  the  NEGF  theory  [10], 


/+oo 

Tr 

-oo 


^ ph 


(E)G>(E)-Z>h(E)G<(E) 


EdE. 


(11) 


The  emission  rate  in  each  mode,  Pg,  can  be  defined  as  Rq  =  Pq/Ti(jOq , 
which  can  be  also  expressed  in  terms  of  absorption  and  emission  processes 
[10],  Aq  and  Eqj 

Rq  =  [(Nq  +  l)Eq-NqAq}.  (12) 

The  calculation  details  of  these  functions  can  be  found  in  [10]  and  [11]. 
Under  stationary  conditions,  equations  (10)  and  (12)  can  be  combined  to 
give  the  explicit  solution 


jy  _  nq(To)Jg  +  Eq 
1  Jq  +  Aq  —  Eq 


(13) 


This  is  a  central  result  of  our  model.  As  an  important  consistence  check,  we 
could  show  that  under  equilibrium  condition  (no  applied  bias),  Aq  —  0  and 
Eq  —  0  and  consequently,  Nq  =  nq(To). 


3.2  Phonon-phonon  decay  rates  into  the  contacts 

The  power  dissipation  from  the  molecule  to  the  contacts  is  described  by 
the  rates  Jq.  These  can  be  computed  from  first  principles  or  treated  as 
parameters.  The  calculation  is  performed  by  setting  up  the  full  Hessian 
and  partitioning  it  in  the  same  way  as  the  Hamiltonian,  namely  into  device 
region  and  contacts.  A  Green’s  function  is  then  built  for  the  vibrational 
eigensystem: 

E  nEj  =  E  (14) 

3  j 
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using  open  boundary  conditions,  as  usual.  In  the  previous  equation  Tiij  is  a 
matrix  element  of  the  Hessian,  e, •  are  the  normal  modes  of  vibration,  Mij  is 
the  mass  matrix  element  and  ujq  is  the  frequency  of  the  mode.  From  (14)  it 
is  possible  to  construct  a  self-energy  for  the  phonons,  whose  imaginary  part 
can  be  used  to  extract  the  phonon  relaxation  lifetime  [11]. 

This  is  essentially  a  Fermi  Golden  Rule  treatment  including  first  order 
one  phonon  to  one  phonon  decay  processes.  It  neglects  a  large  number  of 
other  decay  channels,  for  example  related  to  anharmonic  effects.  Particu¬ 
larly,  high  frequency  modes  do  not  couple  directly  with  the  phonon  band  of 
the  contacts,  meaning  that  only  one  phonon  to  many  phonons  decay  pro¬ 
cesses  are  relevant  in  describing  such  relaxation  mechanisms.  Some  of  these 
effects  can  be  partially  included  by  defining  an  effective  phonon  density  of 
the  contact  modes.  The  formal  Green’s  function  for  the  phonons  is 

G1p^2)  =  Mi -nM-  n^(^2)  -  n^(w2)  (15) 

where  the  self-energies  I\.rL  R{u2)  map  the  infinite  contacts  into  the  finite 
portion  of  the  device.  The  real  part  of  the  self-energies  represent  a  shift 
of  the  frequencies  of  the  normal  modes  of  the  extended  molecule  induced 
by  the  coupling  of  that  region  with  the  bath,  the  imaginary  part  instead 
describes  the  phonon-phonon  decaying  process,  Jq. 


3.3  Definition  of  molecular  temperature 

For  the  purpose  of  our  studies  we  define  an  effective  molecular  temperature, 
Tm,  as  a  direct  measure  of  the  internal  energy  stored  in  the  vibrational 
degrees  of  freedom.  Tm  is  found  by  imposing  the  following  identity, 


U  =  hUqNq 


=  E 


TlUn 


exp( 


flUJ 

kbTm 


- 1 


(16) 


where  U  is  the  total  internal  energy. 

As  shown  in  the  following  section,  as  well  as  in  several  reports  [10,  11,  12], 
this  temperature  parameter  faithfully  maps  the  internal  energy  and  can  be 
considered  as  an  equivalent  measure  of  this  quantity.  We  could  also  think  in 
terms  of  a  nanoprobe  put  in  contact  with  the  molecule  and  able  to  locally 
measure  its  temperature.  Although  this  system  has  not  been  experimentally 
realized  yet,  we  can  think  at  it  as  a  gedanken  experiment.  A  thermometer 
registers  a  temperature  as  a  consequence  of  an  equilibration  process.  Such  an 
equilibration  is  a  consequence  of  an  energy  exchange  between  molecule  and 
thermometer.  The  energy  flux  between  the  two  systems  stops  at  equilibrium, 
that  can  be  expressed  by 

=  EJ<?—  (tp)>  (17) 

q  q  n  0 
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where  Jq  is  the  rate  of  phonon  decays  from  the  molecule  to  the  bath  and 
uq{Tp )  is  the  equilibrium  distribution  of  phonons  in  the  thermometer,  with 
temperature  Tp.  Assuming  the  decay  rates  all  equal,  we  arrive  to  the  con¬ 
clusion  that  Tp  =  Tm. 

3.4  Applications  to  fullerene  on  Cu  and  Si 

The  model  described  above  has  been  applied  to  several  systems.  The  basic 
focus  of  our  investigations  has  been  joule  heating  and  cracking  of  fullerene 
molecules  in  STM  configurations. 

We  have  considered  two  prototype  systems:  Ceo  on  Cu(110)  and  Cgo 
on  Si(100).  These  represent  a  molecule  on  a  metallic  and  a  semiconduct¬ 
ing  surfaces,  and  can  be  useful  to  understand  different  behaviors  due  to 
metallic  bands  and  bandgaps.  The  method  of  choice  for  our  calculation  is 
the  density- functional  tight-binding  method  (DFTB).  The  geometries  of  the 
two  structures,  fully  relaxed  at  the  DFTB  level,  are  shown  in  figure  (6). 
Good  agreement  is  found  with  respect  to  previous  literature. 


Figure  6:  Geometries  of  C6o  relaxed  on  two  different  surfaces,  i)  Si(100):  a) 
top  view  on  the  lowest  fullerene  atoms  close  to  the  silicon  surface,  b)  side  view,  c) 
sketch  of  the  bonding  site,  ii)  Cu(110):  d),  e),  f)  the  same  as  for  silicon. 


The  coupling  between  the  phonons  of  the  substrate  and  the  vibrations 
of  the  molecule  can  be  quantified  by  the  local  density  of  states  (LDOS),  is 
shown  in  figure  (7). 

This  calculations  are  performed  according  to  section  3.2,  using  the  dy¬ 
namical  matrix  (Hessian)  as  produced  by  the  DFTB  code. 

In  the  background  the  surface  density  of  states  for  the  phonons  of  the 
isolated  substrates  is  shown.  The  silicon  substrate  (figure  7(a))  has  a  Debye 
frequency  of  650  cm-1,  while  for  the  much  heavier  copper  atoms  it  is  only 
400  cm-1  (figure  7(b)).  The  broadenings  of  the  molecular  vibrations  are 
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directly  proportional  to  the  decay  rates  into  the  contact  reservoirs,  the  Jq 
coefficient  of  (10).  It  is  easy  to  observe,  that  the  broadening  becomes  very 
narrow  beyond  the  Debye  frequency.  Moreover,  as  expected,  the  broadening 
induced  by  the  silicon  substrate  for  the  modes  below  the  Debye  frequency  is 
substantially  larger  (corresponding  to  a  rate  of  1012  Hz)  than  the  broadening 
for  the  copper  substrate  (1011  Hz).  This  is  a  straightforward  consequence 
of  the  higher  mass  of  copper  compared  to  silicon.  In  order  to  compute  the 


Frequency  (cm  ) 


Frequency  (cm  1 ) 


Figure  7:  The  phonon  LDOS  for  the  two  systems,  a)  Ceo  on  Si(100)  and  b) 
C6o  on  Cu(110).  The  black  line  represents  the  LDOS  with  broadening  induced  by 
coupling  to  the  surface  modes.  In  the  background,  the  shaded  areas  represent  the 
phonon  density  of  states  for  the  isolated  substrates. 


Figure  8:  Calculated  I-V  characteristics,  a)  Ceo  on  Si(100),  b)  Ceo  on  Cu(110). 


current  through  the  system  a  second  contact  has  been  added  to  the  geometry. 
A  copper  STM  tip  was  used  in  both  cases,  In  all  our  simulations  the  substrate 
and  tip  where  assumed  to  be  at  T=0  K.  The  I-V  characteristics  are  plotted 
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in  figure  (8)  for  different  tip-molecule  distances.  The  silicon  substrate  was 
assumed  heavily  p-doped  in  order  to  make  it  conducting,  with  a  Fermi  energy 
at  the  valence  band  edge  with  a  carrier  concentration  of  1019  holes  cm-3.  For 
both  systems  under  investigation  the  tip-molecule  distance  (Z)  and  the  bias 
(V)  control  the  molecular  temperature,  however,  in  the  case  of  the  silicon 
substrate  the  energy  gap  also  plays  an  important  role.  The  temperature  of 
C6o  on  silicon  is  shown  in  figure  (9)  together  with  the  electronic  LDOS  of 
the  molecule.  Due  to  the  energy  gap  the  molecular  states  of  Cgo  couple  only 
weakly  to  the  substrate,  giving  rise  to  a  sharp,  delta-like  shape  of  the  states 
in  the  LDOS.  This  small  broadening  is  reflected  in  the  low  conductivity 
within  the  gap  region. 

Steps  in  molecular  temperature  are  observed  when  molecular  resonances 
enter  the  conduction  window.  However,  due  to  energy  conservation,  the 
electron-phonon  emission  from  (??)  is  related  to  the  Green’s  function  eval¬ 
uated  at  an  energy  shifted  by  uoq.  Hence  the  resonance  is  fully  effective  only 
at  a  bias  beyond  Vres  +  ujmax,  where  Vres  is  the  bias  at  which  the  resonance 
enters  in  conduction  and  ujmax  is  the  largest  frequency  of  the  molecule.  In 
our  fullerene  case,  this  leads  to  a  shift  of  0.2  V  between  Vres  and  the  energy 
position  of  the  temperature  step. 


Figure  9:  C6o  on  Si(100)  (Top)  and  Cu(110)  (Bottom),  a)  temperature  in  the 
junction  for  different  tip-molecule  distances,  b)  internal  energy,  Um,  for  different 
tip-molecule  distances. 


A  similar  behavior  can  be  noticed  for  the  fullerene  on  copper  shown  in 
figure  (9).  Comparing  the  temperature  curves  of  figure  (9)  and  figure  (9), 
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some  important  differences  can  be  noticed.  The  temperature  for  the  copper 
substrate  increases  much  faster  than  the  Si  counterpart. 

One  may  expect  that  the  reason  for  this  is  simply  the  higher  current 
observed  for  the  Cu  substrate,  leading  to  a  larger  joule  heating.  However,  if 
we  carefully  inspect  the  two  curves  we  find  that  this  conclusion  is  deceptive. 
At  1  V  the  C6o  temperature  on  Si  is  700  K  against  1950  K  on  copper  for 
a  tip-molecule  distance  of  2  A.  However,  the  current  is  0.1  fiA  on  silicon 
compared  to  40  /iA  on  copper,  therefore  400  times  smaller.  Based  on  the 
current  argument  we  should  expect  a  much  lower  temperature  for  the  Si 
substrate  (or  a  much  larger  temperature  on  Cu). 

The  reason  for  this  discrepancy  can  be  traced  back  to  electron-hole  (EH) 
pair  excitations,  playing  the  role  of  an  important  damping  mechanism.  In 
metal  this  effect  is  quite  relevant  and  cools  the  molecule  down. 

On  the  other  hand,  on  the  silicon  substrate  the  EH  mechanism  is  inac¬ 
tive  because  of  the  energy  gap  and  the  heating,  despite  a  stronger  bonding 
and  a  larger  phonon  coupling,  is  faster.  Further  important  evidence  for  the 
relevance  of  the  electron-hole  cooling  mechanism  can  be  seen  in  the  temper¬ 
ature  curve  obtained  for  silicon  by  playing  with  the  molecule-tip  distance 
as  shown  in  figure  (9).  Increasing  the  tip-molecule  distance  has  two  oppo¬ 
site  effects.  On  the  one  hand  it  decreases  the  molecular  temperature  due 
to  current  lowering,  on  the  other  hand  it  also  decreases  the  EH  damping 
mechanism,  favoring  heating.  This  explains  why  on  Si,  the  molecular  tem¬ 
perature  tends  to  increase  when  the  tip  moves  from  1.5  A  to  2.0  A.  In  this 
case  the  second  mechanism  prevails  over  the  first. 

The  model  and  the  calculations  presented  above  have  been  also  applied 
to  the  understanding  of  the  cracking  experiment  of  fullerenes  on  Cu  by  the 
approaching  of  an  STM  tip  with  a  fixed  voltage  [13].  In  order  to  understand 
the  observed  I-V  characteristics  before  the  Cgo  irreversibly  deteriorates  we 
have  fitted  the  experimental  data  with  our  model  by  defining  a  critical  tem¬ 
perature  above  which  we  assume  the  fullerene  has  cracked.  This  temperature 
is  found  at  about  1600  K.  From  the  critical  temperature  we  can  extract  a 
critical  tip  distance  and  a  critical  current /volt age.  Figure  10  shows  the  im¬ 
portant  achievements  of  our  theory.  The  critical  current  (I dec  is  plotted  vs 
the  critical  voltage  V dec  and  agrees  well  with  the  experimental  behavior  as 
discussed  in  [13].  An  important  point  is  that  that  our  model  support  the 
explanation  that  the  molecule  can  support  high  currents  when  the  tip  is  in 
contact  because  electron/hole  excitations  provide  a  good  absorption  channel 
keeping  the  molecule  cool. 

3.5  Conclusions 

In  conclusion  of  this  section,  we  are  satisfied  with  the  microscopic  theory 
of  heat  dissipation  that  we  have  implemented.  The  fundamental  parame¬ 
ters  which  have  been  defined  in  the  model  are  the  rate  of  phonons  emitted, 
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Figure  10:  Right,  (a)  Molecular  temperature  due  to  the  flow  of  electronic  current 
as  a  function  of  bias  for  the  indicated  distances,  (b)  The  values  Idee  (triangle, 
green  continuous  line)  and  P dec  (circle,  blue  broken  line)  are  marked  for  a  threshold 
temperature  of  1650  K  of  the  §0  molecule.  The  black  dashed  line  plots  the  position 
of  the  resonances.  The  blue  shaded  area  marks  the  contact  regime,  (c)  Ratio  of 
electron- hole  pairs  in  the  tip  to  total  number  of  electron- hole  pairs  in  the  junction 
as  a  function  of  tip- molecule  distance  (at  0.4  V).  (d)  Temperature  Tm  at  0.4  V  as 
a  function  of  tip-sample  distance  with  (blue)  and  without  (red)  electron-hole  pair 
formation  in  the  tip.  Left.  Experimental  results. 


Eg,  and  the  rate  of  phonon  absorbed,  Ag,  in  each  vibrational  mode  due  to 
electron-phonon  coupling  and  the  rate  of  phonon  dissipated  into  the  heat 
reservoirs,  Jg,  due  to  phonon-phonon  decay  rates.  These  parameters  have 
been  calculated  from  first-principles  in  order  to  account  for  molecular  heat¬ 
ing  as  a  consequence  of  a  current  flux.  Adaptations  could  be  provided  to 
other  mechanisms,  for  which  a  detailed  microscopic  model  should  be  pro¬ 
vided.  For  instance  it  would  be  possible  to  include  sources  of  heat  generated 
by  non-radiative  exciton  recombinations  or  relaxations  of  Frank-Condon  ex¬ 
citations  after  photon  absorption. 

4  Averall  conclusions  and  outlooks 

In  conclusion  we  have  shown  different  heat  transport  approaches,  valid  at 
different  length-scales  ranging  from  macroscopic  domains,  where  the  Fourier 
heat  transport  model  is  appropriate,  to  mesoscale  where  Boltzmann  trans¬ 
port  should  be  applied  down  to  nanoscale  (molecular)  level  where  full  quan¬ 
tum  mechanical  approaches  are  necessary.  In  the  latter  case  we  solve  the 
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quantum  dissipation  within  the  non-equilibrium  Green’s  function  frame¬ 
work.  Different  scales  can  be  coupled  in  a  self-consistent  manner.  The 
bridge  multiscale  approach  must  be  applied  in  this  context,  where  different 
model  domains  must  be  coupled  via  boundary  conditions  on  temperature 
and  energy  fluxes.  As  an  example  of  this  we  have  shown  the  coupling  of 
the  gray  transport  model  and  the  Fourier  model.  As  a  future  development 
of  the  method  we  plan  to  demonstrate  the  coupling  of  the  quantum  dis¬ 
sipation  model  in  molecular  bridges  with  the  gray  model  and  the  Fourier 
model.  At  the  interface  between  quantum  region  and  contact  region  we  will 
impose  the  heat  sources  computed  within  the  nanoscale  approach  and  set 
a  Fourier  or  gray  model  transport  equation  for  a  self-consistent  solution. 
This  approach  will  enable  us  to  evaluate  the  non-equilibrium  temperature 
of  the  nanoscale  contacting  leads,  usually  assumed  in  equilibrium  in  molec¬ 
ular  electronics  calculations.  This  will  give  us  one  of  the  first  examples  of 
a  multiscale/multiphysics  approach  to  the  dissipation  problem  in  nanoscale 
conductors  that  can  cover  length-scales  from  few  nanometers  up  to  microm¬ 
eters. 
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